The quantum phase transition in the sub-ohmic spin-boson model: 
Quantum Monte- Carlo study with a continuous imaginary time cluster algorithm 
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A continuous time cluster algorithm for two-level systems coupled to a dissipative bosonic bath is presented 
and applied to the sub-ohmic spin-Boson model. When the power s of the spectral function /(co) oc co' is smaller 
than 1/2, the critical exponents are found to be classical, mean-field like. Potential sources for the discrepancy 
with recent renormalization group predictions are traced back to the effect of a dangerously irrelevant variable. 
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Quantum-mechanical systems embedded into a dissipative 
environment play an important role in many areas of physics 
111 121] . Among the numerous applications of models that cou- 
ple a small quantum-mechanical system to a bosonic bath are 
noisy quantum dots [3J, decoherence of qubits in quantum 
computations \^, and charge transfer in donor-acceptors sys- 
tems [5']. A major research field are quantum impurity models 
(i.e. a quantum spin embedded in a crystal lattice, for a review 
see ial), where in particular quantum critical points occuiTing 
for instance in the Bose-Fermi Kondo model have been stud- 
ied intensively |7, 8, 9, 10]. 

The paradigmatic model of a two-state system coupled to 
an infinite number of bosonic degrees of freedom is the spin- 
boson model 11,2]. As a function of the strength of the cou- 
pling to its bath it displays a quantum phase transition (QPT) 
at zero temperature between a delocalized phase, which al- 
lows quantum mechanical tunneling between the two states, 
and a localized phase, in which the system ceases to tunnel in 
the low-energy limit and behaves essentially classically. 

While the phase transition is understood in the case of 
ohmic dissipation (s — 1), the sub-ohmic situation (s < 1) 
has been investigated in detail only recently. On general 
grounds, one expects the phase transition to fall into the same 
universality class as that of the classical Ising spin chain 
with long-range interactions \IV\. Indeed, a continuous QPT 
has been found in the spin-boson model for all values of 
< >? < 1 [12], using a generalization of Wilson's numeri- 
cal renormalization group (NRG) technique [6]. However, on 
the basis of these NRG calculations, it was suggested that the 
quantum-to-classical mapping fails for s <\/2 jlSll : There, 
the Ising chain displays a mean-field transition, whereas the 
critical exponents extracted from NRG were non-mean-field- 
like and obeyed hyperscaling. Subsequent NRG calcula- 
tions for the spin-boson [14T and Ising-symmetric Bose-Fermi 
Kondo model \M confirmed this claim. Such a breakdown 
of quantum-to-classical mapping has consequences not only 
for quantum-dissipative systems, but also for Kondo lattice 
models studied within extended dynamical mean-field theory, 
where non-mean-field critical behavior is at the heart of so- 
called local quantum criticality i9[]. 



The purpose of this letter is two-fold: 1) We present a 
novel and accurate quantum Monte-Carlo (QMC) method to 
study the low temperature properties of the sub-ohmic spin- 
boson model, and 2) we determine its critical exponents at 
the quantum phase transition using this method together with 
finite temperature scaling and re-confirm the correctness of 
the quantum-classical mapping for the sub-ohmic bath with 
s<\/2. 

The spin-boson Hamiltonian is defined as 



H^A- 



^ A.; (a,- + a+ ) + ^ CO; a+a; 



(1) 



where a^'^ are Pauli spin- 1/2 operators, a^ , a, are bosonic 
creation and annihilation operators, A the tunnel matrix ele- 
ment, and CO, the oscillator frequencies of the bosonic degrees 
of freedom. The coupling between the the spin o and the bath 
via the A-,- is determined by the spectral function for the bath: 



7(co) = 7t52^?5(K) - CO,) = 27i:a • co^. *co' 



(2) 



for < CO < C0( and 7(co) = otherwise, a represents the cou- 
pling strength to the dissipative bath and co^- is a cut-off fre- 
quency. The parameter s specifies the low-frequency behavior 
of the spectral function: s = \ represents an ohmic bath, and 
,? < 1 a sub-ohmic bath. A system described by ^ and (|2|l 
displays for i < 1 a quantum phase transition (at zero temper- 
ature) at a critical coupling strength a^. In the following we 
determine the critical exponents and herewith the universality 
class of this transition with the help of a continuous time clus- 
ter algorithm that samples stochastically the imaginary time 
path integral for the partition function of the model ([T]i. 

Consider a Hamiltonian for an Ising spin in a transverse 
field of the form 

// = rd'' + G(d',fl,co) , (3) 

where Y is the transverse field strength, Y = A/2 in ([T]), a and 
CO a set of Hermitian operators and parameters, respectively, 
like the Bose operators and coupling constants and frequen- 
cies in the spin-boson model. G is a function of the d~ and (a, 
CO) alone, it is Hermitian but otherwise arbitrary. 



The partition function for this Hamihonian is derived by 
impHcitly performing the Umit of an infinite number of time 
shces in its Suzuki-Trotter representation I15l Il6l llTIl and 
yields the imaginary time path integral 



Z = Tr, 



c,aexp(-Pi/) 

= / Do(t) exp(-54o(T)]) 

P. 



/' 



(4) 
(5) 



where 5m[o(T)]) = -lnTr^exp[/Q iiTG(o(T),a,Co)] and o{l) 
is now a real valued function of the imaginary time T G [0, p], 
denoted as a spin- 1/2 world line. These world lines repre- 
sent realizations of a two-valued Poissonian process that is 
sketched in Fig. la: They are piecewise constant functions 
consisting of consecutive segments of spin-up (a = +1) and 
spin-down (o — —1), where the spin-flips occur at stochas- 
tic times < Ti < Ti,...t„ (n arbitrary) and the interval 
lengths At,- = T,+i — x,- obey a Poissonian statistics P(At) = 
r^' exp(— FAx) with mean value 1/F [ 16]. The path integral 
(|5]l can hence be directly sampled by generating stochastically 
realizations of such world lines and accepting them according 
to their "Boltzmann"-weight exp(— 5o)[a(T)]). More efficient 
sampling procedures like cluster algorithms are based on this 
principle UM- 

For a general transverse Ising model (without coupling to a 
dissipative bath) G(d") represents just the "classical" energy 
E{g^) that is diagonal in the z-representation of the spin- 1/2 
degrees degrees of freedom and 5[a(T)] = /q dzE{o{z)). This 
form holds for an arbitrary number of spins in a transverse 
field, and for arbitrary spin-spin interactions. 

In the case of the spin-boson model ([T]| with the spectral 
function Q the trace over the oscillator degrees of freedom 
yields 101 ■Sm = Ssb with 

SsbHx)] = -I dxj\i'o{i)mi-x')o{x'). (6) 
The kernel imposes long-range interactions in imaginary time: 



Jo 7t 



sinh(fco) 



(7) 



It has the symmetry K{^ — t) = K{i) and the asymptotics 
K{i) oc T"(i+*) for Tc <CT ^ P, where T^ = 271/(0^.. For T < i^ 
the Kernel K{i) is regularized via the frequency cut-off cOt in 
(|2|l and approaches a constant for T -^ 0. 

An efficient way of sampling the path integral is a cluster 
algorithm based on lll6ll . It is generalization of the S wendsen- 
Wang cluster algorithm llSll to continuous time world lines, 
in which not individual spins but the world line segments are 
connected during the cluster-forming procedure, and has to 
incorporate the long-range interactions [19]. It is sketched in 
Fig. lb: Starting from a world line configuration a(T) new 
potential spin-flip sites are introduced according to a Pois- 
sonian statistics, then all segments are pairwise "connected" 
with probability 

pisi,sn) = l-expf-2J^^ d% j^ dx'm%-%'yj . (8) 
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FIG. 1 : (Color online) a Realization of an imaginary time world line 
of a spin- 1/2 in a transverse field, b) Sketch of the continuous time 
cluster update: 1) Starting configuration. 2) Random insertion of new 
potential spin flips (red dots) with Poissonian statistics 3) Connection 
of segments with probabilities given by eq.©. Different colors indi- 
cate the resulting clusters. 4) Each cluster is flipped with probability 
1/2 (the blue one was not flipped). 5) Resulting new imaginary time 
world line. 



where a,b and c,d denote the limits of segment S; and Su, 
respectively. Finally the connected clusters are identified and 
flipped with probability 1/2. All potential spin-flip times that 
do not represent real spin-flips are then removed. We imple- 
mented this algorithm and tested it by comparing results with 
those obtained with conventional Monte-Carlo procedures in 
discrete imaginary time extrapolated to an infinite number of 
time-slices. We analyzed the sampling characteristics of the 
algorithm for the kernel (|7]) with (|2]l over the whole range 
< 5 < 1 and found that on average after 5 updates as sketched 
in Fig. lb the world line configuration are statistically indepen- 
dent from the starting configuration. The data presented below 
represent averages over 10^-10^ cluster updates. 

To study the phase transition in the sub-ohmic spin-boson 
model (s < 0.5) we utilize the finite-^ scaUng forms for ther- 
modynamic observables close to the critical point a = a^ 



(0)r,a = p-^'''^o(p-^'5), 



(9) 



where 5 = (a — ac)/(Xc denotes the distance from the crit- 
ical point, xo and go are the scaling exponent and scaling 
function of the observable O, respectively. The exponent y* 
is 1/v below the upper critical dimension (s > 1/2), v be- 
ing the correlation length exponent, and y,* = 1/y + (1/2 — i) 
above the upper critical dimension (s < 1/2). 11911 . We use the 
dimensionless ratio of moments Q = {m^)^/{m'^), which has 
xq ~ Q and is therefore asymptotically independent of tem- 
perature at 5 = 0, to locate the critical point a^ as shown 
for s — 0.2 in Fig. 2a. This estimate for a^ is then used to 
perform the finite-P scaling analysis for Q, Q — 2(P^' 5) the 
magnetization m — {\m\) — ^^'i'^^m{^^'' 5) and the susceptibil- 
ity X = P(m^) ~ p^-^'A"^%(P^'*5), where y^ is the magnetic ex- 
ponent. The data collapse that one obtains with the mean-field 
values for the exponents yf and y^ 



y; = l/2, yl^ 3/4 



(10) 



is good, as shown Fig.2b-d. 
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FIG. 2: Results for the spin boson model for s = 0.2 and A = 0. 1. a) 
Moment-ratio 2 as a function of the coupling constant a for different 
values of (3. The critical coupling is at a^. = 0.0175 ±0.0002. b-c) 
Finite (3-scaling for the moment-ratio Q, magnetization m and sus- 
ceptibility X (5 = (oc — a.c)/Uc, with ttc from a). The values for the 
critical exponents are y* = 0.5, y\ = 0.75. For large positive values 
of the scaling variable corrections to scaling are stronger. 



At the critical point a = a^ the scaling forms predict % °^ 
T^'^ with X = 2yl — 1 = 1/2, which is clearly confirmed by 
our data displayed in Fig. 2d: % ■ T^'^ coUpaes onto one point 
at 5 = 0. Moreover the scaling forms imply at T = X "^ 
I a — ttc I ^^ with y = {2yl — 1 ) /y,* = 1 , which is demonstrated 
in Fig.3a-c for different values of i < l/2.,andmo= (a — ac)P"' 
for a > a^ with p,„ = —{yl — i)/y't — 1/2^ which is demon- 
strated in Fig. 3d. 

Next we allow for an unbiased fit of the critical exponents to 
our data, including corrections to scaling as in 1 19]. We deter- 
mined y* and yl by finite-P scaling of dQ/da (a = a^ ) ^x p-^v* 
and x(cx = oCc) °= P^^'i^^ The results confirm (fTOl i within the 
error bars for the whole range of s < 1/2 that we studied. 
Only close to 5 = 1/2 the finite-P scaling analysis is impeded 
by the presence of logarithmic corrections at the upper criti- 
cal dimension. Fig.4a-b shows the resulting estimates for the 
exponents 1/v = y* — (1/2 — i) and p,„ = —{yl — l)/.yr* ^s_a 
function of.? in comparison with the NRG predictions of 1 131. 

Although our results for the critical exponents of the sub- 
ohmic bath obtained with our continuous imaginary time algo- 
rithm deviate from the NRG prediction, results for the phase 
diagram match: In Fig.4c our estimates for the critical cou- 
pling ac are compared with those obtained with the NRG 
method [13], they agree very well. 

We confirmed the scenario described here for other values 
of A and (0^, and also for smooth frequency cut-offs as well as 
for other kernels (|7|l, like one that has a regularization in time 
{K{z) = for T < Zc) rather than in frequency. We also found 
that the limit 0)^ — > °° (or T^ -^ 0) exists and is approached 
smoothly and fast, and conclude that, concerning the critical 
exponents, the regularization does not play a significant role. 
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FIG. 3: a-c) Data for the susceptibility X as a function of the distance 
from the critical point 6 > (i.e. in the delocalized phase) for s = 
0.1 (a), s = 0.2 (b) and s = 0.3 for different values of (3 (A and cOc 
as in Fig. 2). For increasing inverse temperature P the data points 
approach the straight line, which is the zero temperature behavior 
X o= 8^'. d) Magnetization m as a function of 5 > (i.e. in the 
localized phase) for p = 2 for different values of i- (multiplied with 
2, 4 and 8 for s = 0.2, s = 0.3 and s = 0.4, respectively, for better 
visibility) The straight lines are guides for the eye proportional to the 
zero temperature behavior (a — ttc)^'' . Shown are only the data 
that are free from finite-P corrections. 



We also implemented a conventional QMC algorithm in 
discrete time (with a finite number of Trotter time slices M) 
and found that for any fixed value of At = ^/M mean-field 
exponents describe perfecdy the scaling at the critical point 
for s < 1/2 (see also il9ll20|] ). Moreover we found that the 
extrapolation M ^ oo of numerical data for Q, m and % ob- 
tained for fixed M reproduces perfectly the results obtained 
with our continuous imaginary time cluster algorithm and that 
the convergence is smooth and fast (with 1/M, as expected). 

Our conclusion therefore is that the quantum-to-classical 
mapping does not fail in the sub-ohmic spin-Boson model. 
The question remains, why the NRG calculation presented in 
I I2II13II yields apparently correct results for quantities like the 
critical coupling, i.e. the phase diagram (see Fig. 4c), but fails 
to predict the correct critical exponents in the case s < 1/2. 

We believe the problem is rooted in a shortcoming of the 
present NRG implementation. As detailed in Ref. |2]J, due 
to the truncation of the bosonic Hilbert space, the NRG - 
while correctly describing the delocalized phase and the criti- 
cal point - it is unable to capture the physics of localized phase 
of the spin-boson model for s < I. Technically, a finite expec- 
tation value (Ot) is accompanied by a mean shift of the bath 
oscillators which diverges in the low-energy limit. Hence, the 
NRG results are expected to be reliable as far as they do not 
involve properties of the localized fixed point. 



The analysis of critical exponents in Ref. [13| now assumed 
that all exponents are properties of the critical fixed point. 
However, this assumption is invalid for the order-parameter 
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FIG. 4: a-b) Numerical estimates of the critical exponents 1/v and 
P„, as a function of s. Triangles: QMC result (obtained from fi- 
nite temperature scaling of the QMC data as described in the text); 
squares: RG results (from |13], compare also |7]); straight lines: 
mean-filed values for s < 1/2. c) Critical coupling strength ttc as a 
function of s diagram for the spin-boson model i ll 121 ) with cut-off fre- 
quency cOf = 1 and tunnel matrix element A = 0.1. Triangles: QMC 
result; squares: NRG results for fixed NRG discretization parameter 
A = 2 (from |12]). Performing the limit A -^ 1 moves the NRG- 
estimates for a^ slighlty downward. 



related exponents p„, and 5„, if the critical fixed point is Gaus- 
sian (like in a (|)^ theory above its upper critical dimension). 
Then, the order parameter amplitude is controlled by a dan- 
gerously irrelevant variable, and p„j and 8m are properties of 
the flow towards the localized fixed point, which in turn is 
not correctly captured by NRG. (Note that 8„j involves the 
non-linear field response at criticality, which is undefined for 
a purely Gaussian theory.) Therefore, the values of P^ and 5„, 
extracted from (present) NRG calculations are unreliable. 

Considering that the NRG calculations nevertheless gave 
well-defined power laws which were moreover consistent with 
hyperscaling, it is worth asking for the underlying reason. We 
conjecture that the artificial Hilbert-space truncation, which 
determines the flow to the "wrong" localized fixed point and 
limits both the field response and the condensate amplitude, 
is equivalent to an operator which is exactly marginal at crit- 
icality in the (|)^ language. Near criticality, this has no con- 
sequences below the upper critical dimension, s > 1/2, as 
the quartic interaction is relevant here, but for s < 1/2 the 
marginal operator instead dominates over the quartic term. 
It is easy to show that an exactly marginal coupling leads to 
yl = {I +-^)/2 —yi with y; — such that hyperscaling is ful- 
filled, while y,* takes its mean-field value - this is what char- 
acterizes the set of NRG critical exponents [13]. (The cor- 
rect result yl = 3/4 for s < 1/2 implies y; = (2^ — l)/4 aris- 
ing from the dangerously irrelevant variable.) The above rea- 
soning is supported by analyzing fermionic impurity models, 
which naturally have the property that a Hilbert space con- 
straint limits the field response. For instance, a resonant-level 
model with power-law bath density of states, which is con- 



trolled by a stable intermediate-coupling fixed point, shows 
hyperscaling for all bath exponents 12211 . 

Finally, the analytical RG argument in Ref.'l3l based on an 
epsilon-expansion for small s, predicted non-mean-field ex- 
ponents obeying hyperscaling for a related reason: While the 
RG equations (9)-(ll) of Ref. 13 are correct, the subsequent 
analysis overlooked the presence of the dangerously irrelevant 
variable, resulting again in the incorrect y^^ = (1 + '^)/2. 

To conclude we have, with the help of an efficient and ac- 
curate continuous time cluster Monte-Carlo algorithm, shown 
that the quantum-to-classical mapping is valid in the case of 
the sub-ohmic spin-boson model. The presence of a danger- 
ously irrelevant variable for s < 1/2 impedes the correct ex- 
traction of the critical exponents with current versions of the 
NRG method - work on its extension to produce reliably the 
necessary determination of magnetic observables in the local- 
ized phase is in progress. 



[1] A. J. Leggett, et al. Rev. Mod. Phys. 59, 1 (1987). 
[2] U. Weiss, Quantum Dissipative Systems (World Scientific, Sin- 
gapore, 1999). 
[3] K. Le Hur, Phys. Rev. Lett. 92, 196804 (2004); M. R. Li, K. Le 

Hur, W. Hofstetter, Phys. Rev. Lett. 95, 086406 (2005). 
[4] M. Thorwart, P Hanggi, Phys. Rev. A 65, 012309 (2001); T. A. 

Costi, R. H. McKenzie, Phys. Rev. A 68, 034301 (2003). 
[5] S. Tornow, N.H. Tong, R. Bulla, Europhys. Lett. 73, 913 (2006). 
[6] R. Bulla, T. A. Costi, T. Pruschke, Rev. Mod. Phys. 80, 395 

(2008). 
[7] M. T. Glossop and K. Ingersent, Phys. Rev. Lett. 95, 067202 

(2005); Phys. Rev. B 75, 104410 (2007). 
[8] M. T. Glossop and K. Ingersent, Phys. Rev. Lett. 99, 227203 

(2007). 
[9] Q. Si, S. Rabello, K. Ingersent, and J. L. Smith, Nature (Lon- 
don) 413, 804 (2001); Phys. Rev. B 68, 115103 (2003). 
[10] L. Zhu, S. Kirchner, Q. Si, A. Georges, Phys. Rev. Lett. 93, 
267201 (2004); J.-X. Zhu, S. Kirchner, R. Bulla, Q. Si, Phys. 
Rev. Lett. 99, 227204 (2007); S. Kirchner, Q. Si, Phys. Rev. 
Lett. 100, 026403 (2008). 
[11] M. Suzuki, Progr Theor. Phys. 56, 1454 (1976). 
[12] R. Bulla, N.H. Tong, and M. Vojta, Phys. Rev. Lett. 91, 170601 

(2003). 
[13] M. Vojta, N.H. Tong, and R. Bulla, Phys. Rev. Lett. 94, 070604 

(2005). 
[14] K. Le Hur, P. Doucet-Beaupre, and W. Hofstetter, Phys. Rev. 

Lett. 99, 126801 (2007). 
[15] C. Pich, A. P. Young, H. Rieger, N. Kawashima, Phys. Rev. 

Lett. 81, 5916(1998). 
[16] H. Rieger, N. Kawashima, Europ. Phys. J. B 9, 233 (1999). 
[17] F. Krzakala, A. Rosso, G. Semerjian, F. Zamponi, 

arXiv:0807.2553 
[18] R.H. Swendsen, J.-S. Wang, Phys. Rev. Lett. 58, 86 (1987). 
[19] E. Luijten, H. W. J. Blote, Phys. Rev. B 56, 8945 (1997). 
[20] M. Troyer (private comminication); S. Kircher and Q. Si (un- 
published). 
[21] R. Bulla, H.-J. Lee, N.H. Tong, and M. Vojta, Phys. Rev. B 71, 

045122(2005). 
[22] L. Fritz and M. Vojta, Phys. Rev. B 70, 214427 (2004). 



